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In finite-size scaling analyses of Monte Carlo simulations of second-order phase transitions one 
often needs an extended temperature/energy range around the critical point. By combining the 
replica-exchange algorithm with cluster updates and an adaptive routine to find the range of in- 
terest, we introduce a new flexible and powerful method for systematic investigations of critical 
phenomena. As a result, we gain two further orders of magnitude in the performance for 2D and 
3D Ising models in comparison with the recently proposed Wang-Landau recursion for cluster 
algorithms based on the multibondic algorithm, which is already a great improvement over the 
standard multicanonical variant. 
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While much attention has been paid in the past to simulations of first-order phase transitions 
and systems with rugged free-energy landscapes in generalized ensembles (umbrella, multicanon- 
ical, Wang -Landau, parallel/simulated tempering sampling) [Jl]], the merits of this non-Boltzmann 
sampling approach also for simulation studies of critical phenomena has been pointed out only 
recently. In Ref. [Q], Berg and one of the authors combined multibondic sampling [J3j] with the 
Wang-Landau recursion [Q] to cover the complete desired "critical" temperature range in a single 
simulation for each lattice size, where the "desired" range derives from a careful finite-size scaling 
(FSS) analysis of all relevant observables. Since the individual reweighting ranges of the involved 
observables may be quite disparate, this scaling analysis is the second important ingredient of the 
method. 

Our new replica-exchange cluster algorithm is a combination of parallel tempering methods [j|] 
with the Swendsen-Wang cluster algorithm For the parallel tempering procedure we use a set 
of A^p replica, where the number of replica depends on the reweighting range that is needed for 
the FSS analysis [0]. To determine this range we perform at the beginning of our simulations a 
short run in a reasonable temperature interval. We choose the number of replica N Tep so that the 
overlap of adjacent histograms is always larger than 25%. This is necessary to ensure that the multi- 
histogram reweighting [Q] works properly. Using the data of this short run as input for the multi- 
histogram reweighting routine, we determine the pseudo-critical points C™ ax = Cl(j8™ x ) of the 
specific heat C(j3) = p 2 V{(e 2 ) - (e) 2 ) and Xl™ of the susceptibility x(P) = PV((m 2 ) - (|m|) 2 ), 
where e = E/V is the energy density and m = M /V the magnetization density. Furthermore, we 
measured the maxima of the slopes of the magnetic cumulants, t/2(j6) = 1 — {m } /3(\m\) and 
£/4(j6) = 1 — (m 4 ) /3(m 2 } 2 , and of the derivatives of the magnetization, d(\m\)/dfi, d(ln\m\)/dfi, 
and d(lnm 2 )/dfi, respectively. We also include the first structure factor (see, e.g., Ref. [^]) in 
our measurement scheme to be directly comparable with the results of Ref. Then we determine 
the j3 values where the observables S = {C,X, • • •} reach the value Sl(P^[ ) = rSf ax , where we 
use the moderate value r = 2/3. This leads to a sequence of values, where (5^ L > /3™ x and 
Asl < fiTl*- I n Fig- 0> we show as an example for such a sequence the reweighted curves for C, 
X, dUz/dfi, and for the two-dimensional (2D) Ising model with linear lattice size L = 8. The 
actual simulation range is then given by the largest interval of the sequence of all fi^i va l ues - m 
our example in Fig. [I], this would lead to an interval [j8^ l>Pcl\- ^ s one can see m ^ s figure, 
the value of L is further away from the critical point then all other ^[ values; therefore, if 
one is not particularly interested in the first structure factor, then the simulation range can be chosen 
narrower. We now use the thus determined interval with the same number of replica for our actual 
measurement run. This interval is usually narrower then the original estimate, hence the overlap of 
the histograms becomes larger and the applicability of the multi-histogram reweighting method is 
assured. 

Let us now illustrate the work flow of our new method for the 2D Ising model. Here we started 
with a reasonable choice of the temperature interval, /3 g ~ = 0.15 and = 0.6, for our smallest 
lattice L = 8. For the large systems we used the measurement interval of the smaller system as 
input interval. Then we used the following general recipe: 

1. choose an input temperature interval 
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Figure 1: Reweighted observables for the 2D Ising model with L = 8. The symbols mark the maximum 
values S™ x and the value S L (j3 5 +/ ~) = nS™ ax with r = 2/3. 



2. choose the number of replica 



3. compute the simulation temperatures for the replica (e.g., equidistant in j8 [10, 11]) 

4. perform several hundred thermalization sweeps 

5. perform a short measurement run 

6. check the overlap of the energy histograms: if the overlap is too small, add 2 replica and goto 
step 3, else go on 

7. use multi-histogram reweighting to determine /3,t L and j8^~ L for all observables S 

8. update the temperature interval according to the largest interval of /3jy and j3^ L , 
i.e., [mins{ps L },max s {p^ L }] 

9. perform several hundred thermalization sweeps 
10. perform the measurement run 

After choosing an input temperature interval and a number of replica for the smallest system, 
our program simulated system sizes from L = 8 up to L = 1024 fully automatically. This shows 
how robust our new method is. Table [I] gives an overview of the automatically determined tem- 
perature intervals, which roughly scale with Lr x l v , where v is the standard critical exponent of the 
correlation length. This scaling can also be used to extrapolate the input interval for larger system 
sizes. In two dimensions, the branch coming from the specific heat has a logarithmically scaling, 
therefore, one could use this knowledge to improve the extrapolation for this special case. We re- 
frain from such modifications to keep the program as generally usable as possible. For comparison 
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Table 1: Simulation range and numbers of replica for the 2D Ising model simulations on L 2 lattices. 
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Figure 2: The temperature interval determined using the specific heat as a function of the system size. The 
horizontal line indicates the critical inverse temperature, the other lines show the exact results calculated 
from the formula of Ferdinand and Fisher [12]. The circles indicate the simulation ranges determined fully 
automatically, cf. Table 1, and the boxes show for completeness the measured values for j3™ ax and fi^ . 



we show in Fig. ^| the calculated temperature interval [P CL ,Pc L ] using the specific heat formula of 
Ferdinand and Fisher [12] and the automatically determined interval of our algorithm. 

To assess the performance of the method, we measured the integrated autocorrelation time 
Ti nt f° r ea °h temperature and system size. We found that the integrated autocorrelation times T; nt 
for the energy, squared magnetization and the first structure factor scale only weakly with the 
system size L. As an example we show Zi nt (E) as a function of L in Fig. |5[ The critical slowing 
down scales oc L z , here we find a dynamical critical exponent z = 0.15(3). When we also take 
the number of replica _/V re p into account and define an effective autocorrelation times T e ff according 
to T e ff = Nrep x Tm t ^ L ZefC , we find a power law with an exponent z e ff = 0.44(2). For T; nt and T e ff 
of m 2 we find slightly smaller values, z = 0.09(3) and Ze,n = 0.37(3), respectively. For the 
dynamical critical exponent is compatible with and for T e ff we find z e ff = 0.29(2). Even the larger 
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Figure 3: Autocorrelation times Tj nt and T e g for the energy of the 2D Ising model, where T e ff = N mp x T^t 
and A^-ep is the number of replica. 



absolute values of the effective autocorrelation times are almost an order of magnitude smaller and 
scale with a much smaller exponent than using the recently proposed multibondic Wang-Landau 
method |@|. 

Due to the fact that in the 2D Ising model the critical exponent a of the specific heat is zero, 
the reweighting range of a single Monte Carlo simulation is oc L _1 / v whereas the range of interest 
scales with L~ r l v (with the range parameter r defined above and v = 1). The number of replica 
needed thus increases with the system size as L( 1_ '^ v , cf. Ref. [13]. InFig.|]we show the numbers 
of replica needed as a function of the system size for various values of r. We also included least 
square fits according to the previously given scaling form and find a reasonable agreement. In 
the 3D Ising model where a«0.11, the reweighting range scales equally to the range of interest 
according to L _1 / v , so that we can use here the same number of replica for all system sizes. In our 
2D simulations only the branch is determined by the scaling of the specific heat. If one omits C 
as a criterion to specify the range of interest in this non-generic case, the numbers of replica can also 
be fixed for all system sizes. As a nice side effect, the dynamical critical exponent for the effective 
autocorrelation times T e ff is even smaller than in the case including C, we find z = z e ff = 0.32(1) for 
the energy and z = z e & = 0.24(1) for m 2 . If the reweighting range is now too narrow to determine 
the critical exponent a directly, one still can use the hyperscaling relation a = 2 — Dv with the 
dimensionality D. 

In Table | we give an overview of the automatically determined temperature intervals for the 
3D Ising model which are similar to the intervals compiled in Table I of Ref. [Q]. By increasing 
the numbers of sweeps in the first short measurement run would lead to a better estimate for the 
temperature interval. We used only about 1 % of our CPU time for this determination, increasing 
this percentage may gain a further improvement of the final results. In Fig. ^| we show the integrated 
autocorrelation times as well as the effective autocorrelation times for the energy of the 3D Ising 
model. Here we find for the dynamical critical exponent z = 0.71(3). In this case T e ff is just a 
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Figure 4: The number of replica needed to cover the range of interest for the specific heat plotted as a 
function of the system size for the 2D Ising model. The straight lines are least square fits according to 
const .lS x ~ r ^l v , with r and v = 1 fixed. 



Table 2: Simulation range and numbers of replica for the 3D Ising model simulations on L 3 lattices. 
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constant shift for all system sizes, due to the fact that the number of replica is independent of 
the system size. We find z = 0.70(3) for the autocorrelation times of m 2 and z = 0.38(4) for 5^. 
In the 3D Ising model the absolute values of the integrated autocorrelation times are almost two 
orders of magnitude smaller and even the effective autocorrelation times are an order of magnitude 
smaller than those reported for the multibondic scheme in Ref. [pp. Since also the dynamical critical 
exponents are smaller, the asymptotic critical slowing down is less pronounced. 

To summarize, we have introduced a very flexible approach for a systematic determination 
and simulation of the critical energy range of interest for second-order phase transitions, which 
one needs to measure critical exponents. The efficiency of the method depends of course on the 
chosen or available update scheme (Metropolis, heat43ath, Glauber, cluster, . . .) in the particular 
case. Since our method is completely general and can be used with any update scheme, it could be 
employed for all simulations in high-energy physics and quantum field theory, statistical physics, 
chemistry and biology where one is interested in critical exponents. 
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Figure 5: T; nt and T e ff for the 3D Ising model, where T e ff = iV rep x Tj nt and N rep = 16 is the number of replica. 
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